The present invention relates to time-of-arrival ranging receivers, and more specifically to global positioning system (GPS) receivers.
A global positioning system (GPS) receiver measures the time of arrival of a satellite ranging code provided by a GPS satellite. Based on the measured time of arrival of the ranging code (which includes orbit parameters and a timestamp indicating when the spacecraft transmitted the ranging code), the receiver determines the time required for the signal to propagate from the satellite to the receiver, and then multiplies this time by the speed of light to obtain a so-called pseudorange; the term pseudorange is used because there are a number of errors involved in the above measurement, the primary error being the error in the time of arrival of the signal because of an offset in the receiver clock from true GPS system time. The true range is just the difference between the time of reception of a specific code phase (of the ranging code) by the GPS receiver and the time of transmission of the specific code phase by the satellite multiplied by the speed of light, i.e.
r=(Tuxe2x88x92Ts)c,xe2x80x83xe2x80x83(1)
where r is the (true) range, Tu is the GPS system time at which a specific code phase is received by the GPS receiver (i.e. the user of the GPS), Ts is the GPS system time at which the specific code phase was transmitted, and c is the speed of light. Because of an offset tu of the receiver clock from GPS system time, what is measured is not the true range r but a pseudorange xcfx81, related to the true range according to,
xcfx81=c[(Tu+tu)xe2x88x92Ts]=r+ctu.xe2x80x83xe2x80x83(2)
The true range r depends on the sought-after user position {right arrow over (u)} (at the time of receipt of the specific code phase) as well as on the position {right arrow over (s)} of the satellite transmitting the specific code phase (i.e. the position of the satellite at the time of transmission of the specific code phase of the ranging code), according to simply,
r=∥{right arrow over (s)}(ttx)xe2x88x92{right arrow over (u)}(trx)∥,xe2x80x83xe2x80x83(3)
where the time of transmission ttx and the time of reception trx are both according to the same clock, the GPS system clock. (FIG. 1 illustrates the notation of equation (3).) Therefore, from equations (2) and (3), the (measured) pseudorange depends on the sought-after user position and receiver clock offset according to,
xcfx81=∥{right arrow over (s)}(ttx)xe2x88x92{right arrow over (u)}(trx)∥+ctu,xe2x80x83xe2x80x83(4)
in which the satellite position {right arrow over (s)}(ttx) is known from the navigation data modulating the ranging code, since the navigation data provide the orbit parameters and a timestamp indicating the time at which the spacecraft transmitted the ranging code, and so indicating the location of the spacecraft along its orbit at the time of transmitting the (specific code phase of the) ranging code.
The problem of determining the offset of the receiver clock is called the time-recovery problem. The prior art teaches solving for the receiver clock offset, as part of the solution for the user position, using various techniques, including linearizing a system of equations for pseudoranges from four or more satellites (using least squares estimation in case of more than four satellites) as well as other techniques such as using a Kalman filter. In a typical solution, at least four satellites are needed because there are four unknowns being solved for all three coordinates of the user position and the receiver clock offset. The approach based on linearizing a system of four pseudorange equations begins with the pseudorange equation, i.e. equation (4), for each of four satellites,   "AutoLeftMatch"                                                                                          ρ                  i                                =                                                      "LeftDoubleBracketingBar"                                                                                                                        s                            ⇀                                                    i                                                ⁢                                                  xe2x80x83                                                ⁢                                                  (                                                      t                            tx                                                    )                                                                    -                                                                        u                          ⇀                                                ⁢                                                  xe2x80x83                                                ⁢                                                  (                                                      t                            rx                                                    )                                                                                      "RightDoubleBracketingBar"                                    +                                      ct                    u                                                                                                                          =                                                                                                                              (                                                                                    x                              i                                                        -                                                          x                              u                                                                                )                                                2                                            +                                                                        (                                                                                    y                              i                                                        -                                                          y                              u                                                                                )                                                2                                            +                                                                        (                                                                                    z                              i                                                        -                                                          z                              u                                                                                )                                                2                                                                              +                                      ct                    u                                                                                                                                            =                                                            f                      i                                        ⁡                                          (                                                                        x                          u                                                ,                                                  y                          u                                                ,                                                  z                          u                                                ,                                                  t                          u                                                                    )                                                                      ,                                                      for                    ⁢                                          xe2x80x83                                        ⁢                    i                                    =                  1                                ,                                  2                  ⁢                                      xe2x80x83                                    ⁢                  …                                ⁢                                  xe2x80x83                                ,                4.                                                                          (          5          )                    
The method the Taylor-series expands each of the equations (5) about an estimated position {right arrow over ({circumflex over (u)})} (related to the sought-after user position {right arrow over (u)} according to {right arrow over (u)}={right arrow over ({circumflex over (u)})}+xcex94{right arrow over (u)}) and an estimated offset {circumflex over (t)}u related to the sought-after receiver offset tu according to tu={circumflex over (t)}u+xcex94tu), and neglects second-order and higher terms in the expansion, leaving,                                           ρ            i                    =                                                    f                i                            ⁢                              xe2x80x83                            ⁢                              (                                                                            u                      ⇀                                        ^                                    ,                                                            t                      ^                                        u                                                  )                                      +                                          ▽                u                            ⁢                                                                    f                    i                                    ⁡                                      (                                                                                            u                          ⇀                                                ^                                            ,                                                                        t                          ^                                                u                                                              )                                                  ·                Δ                            ⁢                              xe2x80x83                            ⁢                              u                ⇀                                      +                                                                                ∂                                          f                      i                                                        ⁢                                      xe2x80x83                                    ⁢                                      (                                                                                            u                          ⇀                                                ^                                            ,                                                                        t                          ^                                                u                                                              )                                                                    ∂                                      t                    u                                                              ⁢              Δ              ⁢                              xe2x80x83                            ⁢                              t                u                                                    ,                              for            ⁢                          xe2x80x83                        ⁢            i                    =          1                ,                  2          ⁢                      xe2x80x83                    ⁢          …                ⁢                  xe2x80x83                ,        4                            (        6        )            
where the notations ∇uƒi({right arrow over ({circumflex over (u)})},{circumflex over (t)}u) and ∂ƒi({right arrow over ({circumflex over (u)})},{circumflex over (t)}u)/∂tu indicate that the function ƒ is to be treated as a function of {right arrow over (u)} and tu for the indicated differentiation, and the result evaluated at {right arrow over ({circumflex over (u)})} and {circumflex over (t)}u. Using equations (5), the equations (6) become,                                           ρ            i                    =                                                    ρ                ^                            i                        -                                                                                                      s                      ⇀                                        i                                    -                                                            u                      ⇀                                        ^                                                                                        r                    ^                                    i                                            ⁢              Δ              ⁢                              xe2x80x83                            ⁢                              u                ⇀                                      +                          c              ⁢                              xe2x80x83                            ⁢              Δ              ⁢                              xe2x80x83                            ⁢                              t                u                                                    ,                              for            ⁢                          xe2x80x83                        ⁢            i                    =          1                ,                  2          ⁢                      xe2x80x83                    ⁢          …                ⁢                  xe2x80x83                ,        4        ,                            (        7        )            
with {circumflex over (ƒ)}i being used in place of ƒi({right arrow over ({circumflex over (u)})},{circumflex over (t)}u) in eq. (6). Using             a      ⇀        i    =                                          s            ⇀                    ^                i            -                        u          ⇀                ^                            r        ^            i      
and xcex94ƒi={circumflex over (ƒ)}ixe2x88x92ƒi, equations (7) can be written as,
xcex94ƒi={right arrow over (a)}ixc2x7xcex94{right arrow over (u)}+cxcex94tu, for i=1, 2 . . . , 4,xe2x80x83xe2x80x83(8)
and by constructing the vectors and matrix,                                                         Δ              ⁢                              xe2x80x83                            ⁢                              ρ                ⇀                                      ≡                                          (                                                                                                    Δ                        ⁢                                                  xe2x80x83                                                ⁢                                                  ρ                          1                                                                                                                                                                        Δ                        ⁢                                                  xe2x80x83                                                ⁢                                                  ρ                          2                                                                                                                                                                        Δ                        ⁢                                                  xe2x80x83                                                ⁢                                                  ρ                          3                                                                                                                                                                        Δ                        ⁢                                                  xe2x80x83                                                ⁢                                                  ρ                          4                                                                                                                    )                            ⁢                              xe2x80x83                            ⁢              H                                =                                                    (                                                                                                    a                        x1                                                                                                            a                        y1                                                                                                            a                        z1                                                                                    1                                                                                                                          a                        x2                                                                                                            a                        y2                                                                                                            a                        z2                                                                                    1                                                                                                                          a                        x3                                                                                                            a                        y3                                                                                                            a                        z3                                                                                    1                                                                                                                          a                        x4                                                                                                            a                        y4                                                                                                            a                        z4                                                                                    1                                                                      )                            ⁢                              xe2x80x83                            ⁢              and              ⁢                              xe2x80x83                            ⁢              Δ              ⁢                              xe2x80x83                            ⁢              x                        =                          (                                                                                          Δ                      ⁢                                              xe2x80x83                                            ⁢                                              x                        u                                                                                                                                                        Δ                      ⁢                                              xe2x80x83                                            ⁢                                              y                        u                                                                                                                                                        Δ                      ⁢                                              xe2x80x83                                            ⁢                                              z                        u                                                                                                                                                                                -                        c                                            ⁢                                              xe2x80x83                                            ⁢                                              t                        u                                                                                                        )                                      ,                            (        9        )            
equations (8) can be written as,
xcex94{right arrow over (xcfx81)}=Hxcex94{right arrow over (x)}xe2x80x83xe2x80x83(10)
which has the solution,
xcex94{right arrow over (x)}=Hxe2x88x921xcex94{right arrow over (xcfx81)}.xe2x80x83xe2x80x83(11)
When navigating by GPS in weak signal conditions such that the navigation component of the navigation signal cannot be decoded, the time when the signal is received is not precisely known. That time, when expressed according to the GPS clock, is called here simply GPS time at arrival and is indicated by the notation xcfx84. It differs from the time at which the signal is transmitted by a satellite (i.e. ttx in equation (5)) by the time of flight tƒ of the signal from the satellite to the GPS receiver. In terms of the GPS time at arrival xcfx84 and the time of flight tƒ, the pseudorange given by equation (2) can be written as,
xcfx81=c└(Tu+tu)xe2x88x92(xcfx84xe2x88x92tƒ)┘=c(Tu+tu)xe2x88x92c(xcfx84xe2x88x92tƒ)=c(tu+tƒ),xe2x80x83xe2x80x83(2xe2x80x2)
since, by definition, Tu=xcfx84.
In normal circumstances (i.e. strong signal conditions), the time at which the ranging code was sent from the spacecraft, xcfx84xe2x88x92tƒ (i.e., in the notation of equation (2), simply Ts), is determined (to within the offset of the satellite clock xcex4t, here taken to be negligible) by decoding the received navigation message. However, when the navigation data (containing the time of transmission/timestamp and the satellite orbit parameters) included in the ranging code sent by the satellite cannot be demodulated (decoded), but the baseband measurements (i.e. aligning the replica code with the transmitted code) can still be made, it is necessary to determine not only the receiver offset tu, but also the time ttx at which the signal was sent by the satellite, which is the GPS time at arrival xcfx84 (the time according to the GPS clock at which the signal was received) minus the unknown time of flight, tƒ.
In the case of signal conditions that are so poor that the navigation component of the ranging code cannot be demodulated, the prior art teaches solving for the user position based on a five-step procedure. The procedure is based on the fact that the pseudorange xcfx81i to the ith spacecraft, which is a measured quantity in the case of strong signal conditions (as opposed to a calculated quantity like the user position), is the sum of two components, and accordingly can be written as,
xcfx81i=xcfx81i(1)+xcfx81i(2)xe2x80x83xe2x80x83(12)
where xcfx81i(1) is obtained from the navigation data included with the ranging code, and is given by                               ρ          i                      (            1            )                          =                  c          ⁢                      ⌈                                          t                f                                                              1                  ⁢                  m                  ⁢                                      xe2x80x83                                    ⁢                  s                                ⁢                                  xe2x80x83                                                      ⌉                                              (        13        )            
where ┌ . . . ┐ indicates the ceiling function (next greatest integer), and xcfx81i(2) is obtained from baseband processing the ranging code (i.e. aligning the replica with the transmitted code), giving the submillisecond part of the pseudorange. The part xcfx81i(1) is here sometimes called the millisecond part, since it gives the pseudorange only to within a millisecond, i.e. to within the speed of light times a millisecond, or nearly 300 km. When the first component xcfx81i(1) cannot be determined (measured) because the navigation data cannot be demodulated, only the second component xcfx81i(2) can be determined (measured). However, with cellular network assistance or some other means of assistance giving us the spacecraft orbit parameters (but not their precise position in the orbits) and the location of a nearby base station, (determined by methods known in the art) the prior art teaches solving for the user position in six steps as follows.
1. Make an estimate of the GPS time at arrival. The estimate can be based on either the local clock or a clock at the reference position.
2. Make estimates of the times of flight for each satellite, using e.g. satellite positions at the estimated GPS time at arrival. If no estimates are otherwise available, use an average of around 70 ms.
3. From the estimates of the GPS time at arrival and times of flight, compute the times of transmission for each satellite and then the satellite positions at those times using orbital parameters for each of the satellites acquired via cellular network assistance.
4. With the satellite positions approximately so determined, use the position of a nearby base station (called here the reference position) to calculate the millisecond part {tilde over (xcfx81)}i(1) of the pseudoranges for each satellite. The full pseudoranges to the receiver (not to base station) are computed by combining the millisecond part {tilde over (xcfx81)}i(1) and the submillisecond part {tilde over (xcfx81)}i(2), where the submillisecond part {tilde over (xcfx81)}i(2) is taken (determined) from the receiver baseband.
5. With the pseudoranges and satellite positions approximately so determined, use the iterative least squares method to determine an estimate of the receiver position. The iterative least squares method gives not only the position, but also an estimate of the a uncertainty in the solution.
6. Repeat steps 1-5 using different estimates of the GPS time at arrival selected from a grid of values around the value provided for example by the cellular network assistance. Use as the user position the value for which the uncertainty is the smallest.
The third step above is quite apparent from the description above. If the time of transmission is known, then the satellite positions are calculated from this time and from orbital parameters using formulas well known in the field of GPS.
In the fifth step, using the spacecraft positions calculated in step three, the user position {right arrow over (u)}(trx) and the receiver clock offset tu are solved using
{tilde over (xcfx81)}i(1)+xcfx81i(2)=∥{right arrow over (s)}i(xcfx84xe2x88x92{tilde over (t)}i(ƒ))xe2x88x92{right arrow over (u)}(trx)∥+ctuxe2x80x83xe2x80x83(15)
At least four such equations are needed, so as to have a system of at least four equations in four unknowns (user position coordinates and receiver clock offset). If we have more than four equations, we can solve the system of equations using least squares, which will give us an uncertainty (error) as part of the solution, and we can iterate on the initial guesses (of first component of the pseudoranges, the times of flight, and the user position) to make the uncertainty smaller, in what is called iterated least squares.
As is evident from the above description, the procedure for solving for user position in case of weak signal conditions is time consuming. It would be advantageous to have a more efficient (faster) way of determining a GPS position-velocity-time (PVT) solution in weak signal conditions, conditions such that the navigation data cannot be demodulated (decoded) but the baseband measurement can still be performed.
Accordingly, the present invention provides a method and corresponding apparatus and system for solving for the position of a satellite positioning system receiver, the receiver having a receiver clock with a receiver clock offset and responsive to ranging codes from some number of satellites, the ranging code from each satellite including a navigation component conveying information about the satellite and having a time of flight from the satellite, the method of use in signal conditions where the ranging codes can be baseband processed so that at least a repeating part of the pseudorange for each satellite can be determined, the method including the steps of: including the system time at arrival as an independent and additional variable of the pseudorange function for each of the ranging codes, an independent variable in addition to the receiver position; providing estimates of the times of flight of the ranging codes and of the repeating part of the pseudoranges for each of the ranging codes; receiving information about the satellite sufficient to determine the orbits of the satellites, although not necessarily the positions of the satellites in their orbits; and solving simultaneously for the receiver position and the system time at arrival.
In a further aspect of the invention, the method also includes the receiver clock offset as an independent and additional variable. In some such applications, the method also includes the step of performing a Taylor-series expansion of the pseudorange function about an initial guess or estimate of each of the independent variables for as many of the satellites as there are independent variables thereby providing a system of pseudorange equations. The pseudorange equations are sometimes solved using a determinate single-point solution that is either an iterative least squares solution or an exact algebraic solution, depending on the number of pseudorange equations, and sometimes solved using a Kalman filter.
In another further aspect of the invention, the information about each satellite includes the orbit parameters for each satellite.
In yet another further aspect of the invention, the satellite positioning system is the global satellite positioning system (GPS) and the system time at arrival is the GPS time at arrival.
In yet still another further aspect of the invention, the repeating part of the pseudorange is the submillisecond part of the pseudorange.
In yet still even another further aspect of the invention, the method is of use in signal conditions where the ranging codes can be baseband processed providing at least the submillisecond part of the pseudorange.
Compared to the prior art procedure explained above, the present invention solves for the user position using the following steps, of which the first four are the same as in the prior art.
1. Make an estimate of the GPS time at arrival. The estimate can be based on either the local clock or a clock at the reference position.
2. Make estimates of the times of flight for each satellite, using e.g. satellite positions at the estimated GPS time at arrival. If no estimates are otherwise available, use an average of around 70 ms.
3. From the estimates of the GPS time at arrival and the times of flight, compute the times of transmission for each satellite and then the satellite positions at those times using orbital parameters for each of the satellites acquired via cellular network assistance.
4. With the satellite positions approximately so determined, use the position of a nearby base station (called here the reference position) to calculate the millisecond part {tilde over (xcfx81)}i(1) of the pseudoranges for each satellite. The full pseudoranges to the receiver (not to base station) are computed by combining the millisecond part {tilde over (xcfx81)}i(1) and the submillisecond part {tilde over (xcfx81)}i(2) where the submillisecond part {tilde over (xcfx81)}i(2) is taken (determined) from the receiver baseband.
5. With the pseudoranges and satellite positions known, instead of using the iterative least squares method as in the prior art, include in the pseudorange equations the unknown GPS time at arrival xcfx84 as an additional, independent variable. Then use the iterative least squares technique for the modified equations, the iteration here being different than in the prior art. In the conventional solution of the prior art, the satellite positions are fixed while iterations refine the estimate of the user position. Here the satellite positions and the user position are refined at the same time. The satellite positions are refined by iterating on the GPS time at arrival xcfx84, which amounts to adjusting the positions of the satellites in their orbits, with the orbits themselves known via cellular network assistance or some other source. (New estimates of the GPS time at arrival and previous estimates of the times of flight are used to compute new times of transmission. The new times of transmission are used to compute new satellite positions. Then the times of flight and pseudoranges are updated.) Each iteration gives a new estimate of the user position, the GPS time at arrival, and the receiver clock offset. The iterations here converge, i.e. the solutions stop changing, and so the iteration here is fundamentally different than that of the prior art.
In essence, the invention collapses the last steps of the solution of the prior art into one step by introducing the GPS time at arrival xcfx84 (the time at which the user receives the ranging code) as a fifth unknown in the pseudorange equation,
xe2x80x83{tilde over (xcfx81)}i=∥{right arrow over (s)}i(xcfx84xe2x88x92{tilde over (t)}i(ƒ))xe2x88x92{right arrow over (u)}(trx)∥+ctuxe2x80x83xe2x80x83(16)
where the pseudoranges {tilde over (xcfx81)}i and the times of flight {tilde over (t)}i(ƒ) are (approximately) known by some means. Thus, the invention provides five unknowns (xcfx84, ux, uy, uz, tu) which are solvable by at least five pseudorange equations (i.e. equation (16) for five satellites). In case of five such pseudorange equations, the invention provides either a determinate solution or a Kalman filter solution. In case of more than five pseudorange equations, the invention provides a solution using either the method of least squares or a Kalman filter.
The method provided by the invention is faster than methods known in the art (from five to ten times faster), and it produces clock times (receiver offsets) that are as accurate as at least some of the prior art methods (in particular, the so-called SSE method). In the method using a filter such as a Kalman filter, the time error can be determined more precisely and the extra computation load is minimized when the time error solution is integrated into the Kalman filter.